clear; clc;
% 固定 A 扫描 f，绘制 bifurcation diagram 和 LE vs f

% 参数设定
a = 0.7; 
b = 0.8; 
c = 0.1; 
alpha = 0.25; 
beta = 0.01; 
delta = 0.1;
% 初始状态 [x; y; z]
y0 = [0.2; 0; 0.01];

T_total = 1200;   % 总模拟时间
T_trans = 400;    % 舍去瞬态时间
dt = 0.01;        % 用于 LE 估计的重新归一化间隔
h = 0.01;         % RK4 积分步长

% 固定 A = 0.35，扫描频率 f
A_fixed = 0.35;
f_range = linspace(0.1, 1.2, 600);

% 用于存储 bifurcation 数据 (f, x_sample)
bifurcationF = [];
% 用于存储 LE 值
LE_f = zeros(size(f_range));

for i = 1:length(f_range)
    f_val = f_range(i);
    params = struct('a', a, 'b', b, 'c', c, 'alpha', alpha, ...
                    'beta', beta, 'delta', delta, 'A', A_fixed, 'f', f_val);
                
    % 数值积分，使用自写的 RK4 积分器
    tspan = [0, T_total];
    [t, Y] = rk4(@(t,y) hybridNeuron(t, y, params), tspan, y0, h);
    
    % 舍去瞬态数据
    idx = t > T_trans;
    t_post = t(idx);
    Y_post = Y(idx,:);
    
    % 根据外部激励周期采样：每个周期采集 x 的峰值
    T_period = 2*pi / f_val; 
    x_samples = [];
    % 注意：ceil(T_trans/T_period) 到 floor(T_total/T_period) 为采样周期数
    for k = ceil(T_trans/T_period):floor(T_total/T_period)
        t_sample = k * T_period;
        [~, idx_sample] = min(abs(t - t_sample));
        x_samples(end+1) = Y(idx_sample,1);
    end
    
    % 累加当前 f 值对应的分岔数据
    bifurcationF = [bifurcationF; repmat(f_val, length(x_samples), 1), x_samples'];
    
    % 计算 LE（使用自写的 RK4 积分器）
    LE_f(i) = LE(params, y0, T_total, dt, h);
end

% 左图：分岔图（x_peak vs f）
subplot(1,2,1);
plot(bifurcationF(:,1), bifurcationF(:,2), 'r.', 'MarkerSize', 5);
xlabel('Frequency f');
ylabel('x_{peak}');
title('Bifurcation diagram (A = 0.35)');

% 右图：LE vs f
subplot(1,2,2);
plot(f_range, LE_f, 'g-', 'LineWidth', 1.5); hold on;
yline(0, 'r--');  % 绘制 LE=0 的参考线
xlabel('Frequency f');
ylabel('LE');
title('LE vs f (A = 0.35)');
